% Aproximace kaskad: Glauber-Sudarshanova funkce aproximovana gaussovskym rozdelenim amplitudy,
% exp(-(r-r0)^2/(2sigma^2)), hledaji se parametry r0 a sigma aby byly splnene hodnoty mean(n) a mean(n^2)
% numericka integrace v 2D

clear
global Nmean N2mean nbar

% Definice funkci
%%%%%%%%%%%%%%%%%%%%%%
% pro nalezeni r0 a sigma z nmean a n2mean
function y=OdchylPar(x)
 global Nmean N2mean
 r0p=x(1);
 sigp=x(2);
 u=r0p/sigp/sqrt(2);
 nmeanp=sigp^2/2*(sqrt(pi)*u*(2*u^2+3)*(erf(u)+1)+2*exp(-u^2)*(u^2+1))/(exp(-u^2)+sqrt(pi)*u*(erf(u)+1));
 citatel=sqrt(pi)*u*(4*u^2*(u^2+5)+15)*(erf(u)+1)+2*exp(-u^2)*(2*u^4+9*u^2+4);
 jmenovatel=sqrt(pi)*u*(erf(u)+1)+exp(-u^2);
 n2meanp=sigp^4/4*citatel/jmenovatel +nmeanp;
 y=(nmeanp-Nmean)^2+(n2meanp-N2mean)^2;
endfunction

% pro nalezeni optimalniho chi a alpha
function y=Dolik(x)
 global nbar
 alpha=x(1);
 chi=x(2);
 G=cos(alpha)^2*sin(alpha)^2*sin(chi)/nbar^2/((1/nbar+2*sin(alpha)^2*sin(chi/2)^2)^2+sin(alpha)^4*sin(chi)^2)^2;
 % koherentni stavy
 %G=nbar*cos(alpha).^2.*exp(-4*sin(alpha).^2*nbar.*sin(chi/2).^2).*besselj(1,2*sin(alpha).^2*nbar.*sin(chi));
 y=-cos(alpha)^2*nbar-G;
endfunction


% Pocatecni hodnoty mean n a mean n^2
Nmean0=50
%N2mean0=Nmean0^2+Nmean0 % Coherent
N2mean0=2*Nmean0^2+Nmean0 % Thermal
%N2mean0=14030

Nmean=Nmean0;
N2mean=N2mean0;

relatfluct0=sqrt(N2mean0-Nmean0^2)/Nmean0

dx=0.5;
xmax=55;
x=0:dx:xmax;
[X,Y]=meshgrid(x,x);


PocetKaskad=10
alphas=zeros(1,PocetKaskad);
chis=zeros(1,PocetKaskad);
Nmeans=zeros(1,PocetKaskad);
N2means=zeros(1,PocetKaskad);
nmeans=zeros(1,PocetKaskad);
n2means=zeros(1,PocetKaskad);
relatflucts=zeros(1,PocetKaskad);
pomnormals=zeros(1,PocetKaskad);
pomrs=zeros(1,PocetKaskad);
pomsigmas=zeros(1,PocetKaskad);
rozdilyNs=zeros(1,PocetKaskad);
rozdilyN2s=zeros(1,PocetKaskad);
chybyoptimalizace=zeros(1,PocetKaskad);

arpg0=[0.25;0.2];
arpg=arpg0;
pomodhad=[sqrt(Nmean),sqrt(Nmean)];


for n=1:PocetKaskad

 vals=fminunc('OdchylPar',pomodhad);
 pomodhad=vals;
 vals=fminunc('OdchylPar',pomodhad);
 pomodhad=vals;
 vals=fminunc('OdchylPar',pomodhad);
 pomodhad=vals;
 vals=fminunc('OdchylPar',pomodhad);
 pomodhad=vals;
 
 r0=vals(1);
 sigma=vals(2);
 pomrs(n)=r0;
 pomsigmas(n)=sigma;
 chybyoptimalizace(n)=OdchylPar(vals);

 u=r0/sqrt(2)/sigma;
%A=1/2/pi/sigma^2/(exp(-r0^2/2/sigma^2)+sqrt(pi/2)*r0/sigma*(erf(r0/sqrt(2)/sigma)+1))
 nmean=sigma^2/2*(sqrt(pi)*u*(2*u^2+3)*(erf(u)+1)+2*exp(-u^2)*(u^2+1))/(exp(-u^2)+sqrt(pi)*u*(erf(u)+1));
 citatel=sqrt(pi)*u*(4*u^2*(u^2+5)+15)*(erf(u)+1)+2*exp(-u^2)*(2*u^4+9*u^2+4);
 jmenovatel=sqrt(pi)*u*(erf(u)+1)+exp(-u^2);
 n2mean=sigma^4/4*citatel/jmenovatel +nmean;
 rozdilyNs(n)=Nmean-nmean;
 rozdilyN2s(n)=N2mean-n2mean;

 % Optimalizace chi a alpha (pro termalni stav)
 nbar=nmean;
 vals=fminunc('Dolik',arpg);
 arpg=vals;
 alpha=vals(1);
 chi=vals(2);
 %%% POZOR - neoptimalizovane hodnoty
 alpha=0.19;
 chi=0.2;
 alphas(n)=alpha;
 chis(n)=chi;
 c=cos(alpha);
 s=sin(alpha);

 A=1/2/pi/sigma^2/(exp(-r0^2/2/sigma^2)+sqrt(pi/2)*r0/sigma*(erf(r0/sqrt(2)/sigma)+1));
 g0=X.*Y.*exp(-((X-r0).^2+(Y-r0).^2)/2/sigma^2);
 g1=X.*Y.*(X.^2+Y.^2).*exp(-((X-r0).^2+(Y-r0).^2)/2/sigma^2);
 g2=X.^2.*Y.^2.*exp(-((X-r0).^2+(Y-r0).^2)/2/sigma^2).*exp(-(s^2*(X.^2+Y.^2)*sin(chi/2)^2)).*besselj(1,s^2*X.*Y*sin(chi));
 gtot=g1+2*g2;

 pomnormal=4*pi^2*A^2*sum(sum(g0))*dx^2;
 pomnormals(n)=pomnormal;

 nmeannew=pi^2*A^2*c^2*sum(sum(gtot))*dx^2;

 f1=(X.^4+Y.^4+4*X.^2.*Y.^2).*X.*Y.*exp(-((X-r0).^2+(Y-r0).^2)/2/sigma^2);
 f2=4*X.^2.*Y.^2.*(X.^2+Y.^2).*exp(-((X-r0).^2+(Y-r0).^2)/2/sigma^2).*exp(-(s^2*(X.^2+Y.^2)*sin(chi/2)^2)).*besselj(1,s^2*X.*Y*sin(chi));
 f3=2*X.^3.*Y.^3.*exp(-((X-r0).^2+(Y-r0).^2)/2/sigma^2).*exp(-(s^2*(X.^2+Y.^2)*sin(chi)^2)).*besselj(2,s^2*X.*Y*sin(2*chi));
 ftot=f1+f2+f3;
 n2meannew=c^4*pi^2*A^2/4*sum(sum(ftot))*dx^2+nmeannew;

 relatflucnew=sqrt(n2meannew-nmeannew^2)/nmeannew;
 
 Nmean=nmeannew;
 Nmeans(n)=Nmean;
 N2mean=n2meannew;
 N2means(n)=N2mean;
 relatflucts(n)=relatflucnew;
end

velik=20;

figure(2)
plot([0:PocetKaskad],[Nmean0,Nmeans],'.')
set(gca,'FontSize',velik)


figure(9)
plot([0:PocetKaskad],[relatfluct0,relatflucts],'.')
set(gca,'FontSize',velik)

figure(3)
plot([1:PocetKaskad],alphas,'.')
%plot([1:PocetKaskad],N2means,'.')
set(gca,'FontSize',velik)

figure(4)
plot([1:PocetKaskad],chis,'.')
set(gca,'FontSize',velik)

figure(5)
plot([1:PocetKaskad],pomnormals,'.')
set(gca,'FontSize',velik)

figure(6)
plot([1:PocetKaskad],[pomrs;pomsigmas],'.')

figure(7)
plot([1:PocetKaskad],[rozdilyNs;rozdilyN2s],'.')

figure(8)
plot([1:PocetKaskad],chybyoptimalizace,'.')
set(gca,'FontSize',velik)